Holstein polaron in two and three dimensions by quantum Monte Carlo 
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A recently developed quantum Monte Carlo approach to the Holstein model with one electron 
[PRB 69, 024301 (2004)] is extended to two and three dimensional lattices. A moderate sign 
problem occurs, which is found to diminish with increasing system size in all dimensions, and not 
to affect simulations significantly. We present an extensive study of the influence of temperature, 
system size, dimensionality and model parameters on the small-polaron cross over. Results are 
extrapolated to remove the error due to the Trotter discretization, which significantly improves the 
accuracy. Comparison with existing work and other quantum Monte Carlo methods is made. The 
method can be extended to the many-electron case. 

PACS numbers: 63.20.Kr 71.27.+a 71.38.-k 02.70.Ss 



INTRODUCTION 



In recent years, a lot of experimental evidence has 
been given for the importance of electron-phonon interac- 
tion in strongly correlated systems such as the cuprates^ 
or the manganitesi^ Although considerable theoretical 
progress has been made in understanding and describ- 
ing many details of the physics of these compounds, 
the quantum nature of the phonons has often been ne- 
glected in actual calculations However, a quantum- 
mechanical treatment of the lattice degrees of freedom 
has been shown to give significantly different results in 
certain parameter regimes;— 

Due to the complexity of the models for the above- 
mentioned classes of materials, numerical methods have 
been used extensively. One very powerful approach is the 
quantum Monte Carlo (QMC) method, which allows for 
simulations on relatively large lattices and gives quasiex- 
act results (i.e., exact apart from statistical errors which 
can, in principle, be made arbitrarily small) also at finite 
temperature. The latter point represents a major advan- 
tage over, e.g., the density matrix renormalization group 
(DMRG) or variational methods — which are restricted 
to the calculation of ground-state properties — since fas- 
cinating phenomena such as high-temperature supercon- 
ductivity and colossal magnetoresistance can be investi- 
gated. Exact diagonalization (ED) can be performed at 
finite temperatures (see, e.g., Ref. but for the case of 
electron-phonon systems such calculations would be re- 
stricted to very small systems and unrealistic parameters. 
Consequently, all existing work based on ED is for zero 
temperature only. Additionally, for coupled electron- 
phonon systems, the infinite-dimensional Hilbert space 
associated with the boson degrees of freedom represents 
a substantial difficulty for ED and DMRG, in contrast to 
QMC. Nevertheless, QMC methods are often limited by 
(1) the minus-sign problem, which restricts simulations 
to high temperatures and/or small systems, (2) the fact 
that the calculation of dynamical properties, such as the 
one-electron Green function, requires analytic continua- 
tion to the real-time axis which is an ill-posed problem. 



and (3) by strong autocorrelations and large statistical 
errors. 

In a recent paper, ^ from here on referred to as I, 
we have proposed a new QMC approach to the Hol- 
stein model, which is based on the canonical Lang-Firsov 
transformation^ and a principal component representa- 
tion of the phonon degrees of freedom. The resulting 
algorithm has been shown to completely avoid the prob- 
lem of autocorrelations, strongly reduce statistical errors 
and to allow for very efficient simulations i& The work 
here and in I is restricted to the case of a single electron 
which is often called the Holstein polaron prohlern!^ . Ob- 
viously, in connection with materials such as the cuprates 
or manganites, we are interested in studying the many- 
electron system. However, the one-electron case repre- 
sents an important first step towards an understanding 
of more general situations. Moreover, due to the coupling 
to the phonons, even the model with one electron consti- 
tutes a complex many-body problem which has received 
a lot of attention in the past. An important feature of 
our approach is the fact that it can be generalized to the 
more realistic many-electron case encountered in transi- 
tion metal oxides. 

In this paper, we extend the QMC method proposed 
in I to square (cubic) lattices in two (three) dimensions. 
Although the one-dimensional case considered in I is of 
great theoretical interest, real strongly correlated materi- 
als usually require two or three dimensional models. Con- 
sequently, before turning to the more challenging many- 
electron case, it is crucial to examine the influence of 
dimensionality on our algorithm. To this end, we inves- 
tigate the dependence of the minus-sign problem, which 
has been found to appear in simulations of the Lang- 
Firsov transformed model, on the dimensionality D of the 
lattice, as well as on the parameters of the model. We 
then use our method to study in detail the well-known 
small-polaron cross over, which occurs with increasing 
electron-phonon coupling strength, and its dependence 
on phonon frequency, temperature, system size and di- 
mensionality. We remove the Trotter error, which results 
from the discretization of imaginary time, by extrapolat- 
ing our results to the limit of continuous time, thereby 
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significantly improving the accuracy of the data. 

The paper is organized as follows. In Sec. ^ we review 
the physics of the Holstein model with one electron, and 
in Sec. lIIII we briefly discuss the QMC method introduced 
in I. Results arc presented in Sec. II VI and Sec.lVlcontains 
our conclusions. 

II. THE HOLSTEIN MODEL 

The Holstein Hamiltonian^ with dimensionless phonon 
variables reads^ 

^ K + P + I, (1) 

where we have introduced the abbreviations K , P and / 
for the kinetic, potential and interaction terms, respec- 
tively. Spin indices have been suppressed in the notation, 
owing to the spin symmetry of the one-electron problem. 
In Eq. J^l, c\ (ci) creates (annihilates) a spinless electron 
at lattice site z, Xi and pi denote the displacement and 
momentum operator of a harmonic oscillator at site i, 
and fii ~ c\c^. The coupling term / describes the local 
interaction of the single electron considered here with dis- 
persionless Einstein phonons. In the first term, the sym- 
bol {ij) denotes a summation over all hopping processes 
i ^ j and j ^ i between neighboring lattice sites 
The parameters of the model are the hopping integral t, 
the phonon energy lu {h = 1), and the electron-phonon 
coupling constant a. As in I, we define the dimension- 
less coupling constant A = a'^/{uiW), where W = 4tD is 
the bare bandwidth in D dimensions. We shall also use 
the dimensionless phonon frequency lu = uj/t^ also called 
the ^^adiabatic ratid\ and express all energies in units 
of t. Consequently, the model can be described by two 
independent parameters, u) and A. While the bulk of re- 
sults in this paper has been calculated assuming periodic 
boundary conditions in real space, we shall also discuss 
the effect of the boundary conditions on the sign problem 
in Sec. HvH 

Since the literature on the Holstein polaron problem 
is vast, we restrict our discussion to work in more than 
one dimension. A brief review of the one-dimensional 
model can be found in I. Moreover, we will focus on re- 
cent progress in the field, and on numerical methods. 
A comprehensive overview of earlier analytical work can 
be found, e.g., in the books of Alexandrov and Mot^ 
and Mahan.^° The Holstein polaron in D > 1 has been 
studied using a large variety of numerical techniques. In 
contrast to many perturbative approaches the lat- 
ter can also accurately describe the physically most in- 
teresting regime of small but finite phonon frequency 
(0 < [D < 1), and intermediate electron-phonon cou- 
pling (A « 1). Much information about the Holstein 
polaron has been obtained using QMC. De Raedt and 
LagendijkiiiiSii^ applied Feynman's path-integral tech- 
nique to the lattice problem of Eq. . In this approach. 



the phonon degrees of freedom are integrated out analyt- 
ically, and one is left with a single-particle system with 
a retarded self-interaction, which can be simulated using 
the Monte Carlo method. The only approximation con- 
sists of discretizing the imaginary time using the Suzuki- 
Trotter decomposition^ As the phonons have been com- 
pletely eliminated from the problem, very efficient sim- 
ulations on large lattices even in higher dimensions are 
possible. Moreover, the method is not restricted to the 
Holstein Hamiltonian (Q. It can be extended to include 
long-range electron-phonon coupling as well as disper- 
sive phonons.^^ Kornilovitcbi^ later employed the same 
method to study the whole range of values of the adi- 
abatic ratio w, and extrapolated the results to remove 
the error due to the Trotter approximation. He also 
pointed out the connection of the algorithm to world- 
line methodsi^ Extending this work, Kornilovitch devel- 
oped a QMC approach which is formulated in continuous 
imaginary time, and which is capable of directly measur- 
ing the polaron band dispersion E{k) and the density of 
states in one to three dimensions, by sampling world lines 
with open boundary conditions in imaginary time^^ii^ 
Although the method gives very accurate results and can 
also be applied to more general models with, e.g., long- 
range interaction, it is limited to the regime of intermedi- 
ate and strong electron-phonon coupling as well as (D > 1 
by a minus-sign problem. ^^'^^ Consequently, it is not as 
universally applicable as the algorithm of de Raedt and 
LagendijkiiiiiSiiSiii As the above discussion reveals, the 
QMC methods of Refs. H 12 13 14 15- 16 are very well 
suited to study the one-electron problem. Moreover, de 
Raedt and Lagendijk also extended their approach to the 
bipolaron problem of two electrons with opposite spins^ 
However, as pointed out by Kornilovitch, » these are all 
world-line methods. Consequently, despite the possibility 
of integrating out the phonon degrees of freedom even in 
the many-electron case, they face a serious sign problem 
in more than one dimension, for two or more fermions 
of the same spin, and can therefore not be used to study 
many-particle systems. In the context of superconductiv- 
ity, the Holstein and the Holstein-Hubbard model with 
many electrons have been investigatedi2ii2iS2i2iiS^ using 
the grand-canonical determinant method of Blankenbe- 
cler et ali^ (see also I). However, as discussed in I, due to 
strong autocorrelations, these simulations were restricted 
to rather large phonon frequencies uj > 1, while, e.g., 
the cuprates and the manganites fall into the adiabatic 
regime a) <^ 1. This is exactly the point where our new 
approach enters. As we shall see below, it completely 
avoids the problem of autocorrelations. 

Apart from QMC, several other methods have been 
successfully applied to the Holstein polaron. This in- 
cludes ED in combination with a truncation of the 
phonon Hilbert spacejSiSiiS^ finite-cluster strong cou- 
pling perturbation theory (SCPT)f2& cluster perturba- 
tion theory (CPT)^ DMRG^ and a variational diag- 
onalization method«2S As pointed out before, standard 
ED is restricted to rather small numbers of lattice sites 



3 



and/or phonon states, while DMRG^^ and the variational 
method of Ref.^^are applicable over a wider range of val- 
ues of phonon frequency and electron-phonon coupling, 
and are much less influenced by finite-size effects. The 
same is true for SCPT and CPT, which exactly diagonal- 
ize small clusters — for which enough phonon states can 
be included in the calculation — and extrapolate the re- 
sults to the thermodynamic limit by treating the rest of 
the system as a perturbationiSSi2i Nevertheless, as men- 
tioned before, all this work was limited to zero temper- 
ature. Finally, the Holstein polaron has been investi- 
gated recently using weak- and strong-coupling pertur- 
bation theory)22i^ and a variational wave function which 
is a superposition of Bloch states for the weak- and the 
strong-coupling regime 

From all this work, many properties of the Holstein 
polaron arc well understood. Similar to one dimensiouf^ 
there is a cross over from a quasiparticlc with slightly 
increased effective mass to a heavy small polaron as the 
electron-phonon coupling strength increases. As pointed 
out by Fehske et alw^ and Capone et al.^ the two con- 
ditions for the existence of a small polaron are A > 1 and 
AD/(Ii > 1. We will first discuss the weak-coupHng state. 
Concerning its nature in higher dimensions, there exist 
two different views. From calculations based on the adi- 
abatic approximation, i.e., taking the limit — > 0, one 
expects a qualitatively different behavior in one dimen- 
sion compared to D > 1, also for tD > 0. Wellein et alS^ 
distinguish between the adiabatic {Q < 1) and the nona- 
diabatic {Co > 1) regime. In the adiabatic case, and for 
D > 1, the electron is expected to remain quasifree with 
an almost unchanged effective mass and kinetic energy 
for A < 1, corresponding to a quasiparticlc with infinite 
radius. This contrasts strongly with the one-dimensional 
case, in which the electron is always self-trapped by the 
surrounding lattice distortion and forms a polaron with 
finite radius for any A > 0. On the other hand, in the 
nonadiabatic regime (D > 1, the behavior is very similar in 
all dimensions, and a very gradual decrease of the kinetic 
energy (or increase in effective mass) is observed as the 
coupling strength increases. Due to the large energy of 
the phonon excitations, only the zero-phonon state con- 
tributes significantly. Moreover, the phonons are fast and 
react almost instantly to the motion of the electron. Con- 
sequently, a lattice distortion only persists in the imme- 
diate vicinity of the electron, and this rather small quasi- 
particlc is sometimes called a ^^nonadiabatic Lang-Firsov 
polaron"^ Romero et al^ take on a slightly different 
viewpoint. They argue that the large polaron state is 
essentially the same in any dimension, and that the only 
effect of increasing the dimension of the system is given 
by the observed sharpening of the transition to a small 
polaron in D > 1 compared to ID. For A larger than 
a critical value Ac, determined by the aforementioned 
conditions A > 1 and AD/u) > 1, where the cross over 
occurs, both views agree on the existence of a a small- 
sized so-called ^'Holstein polaron" ?^ The latter is a heavy 
quasiparticlc with a strongly reduced mobility. The cross 



over at Ac is very sharp, especially for (D ^ 1, but it does 
not represent a real phase transitioni^ Even though the 
electron is trapped in the potential well originating from 
the response of the lattice to its motion, the ground state 
is still Bloch-like. For simplicity, in the remainder of this 
paper, we shall always refer to the weak-coupling state as 
a large polaron, either with finite or infinite radius, de- 
pending on which of the abovementioned viewpoints one 
holds. In our opinion, a definite decision about which 
of the two alternatives is correct cannot be made using 
QMC, which is restricted to finite clusters and, more im- 
portantly, finite temperatures. As a consequence, there 
will always be a contribution from excited states, making 
it difficult to reveal the true nature of the ground state 
in the weak-coupling regime. 



III. QUANTUM MONTE CARLO METHOD 

Since the extension of the QMC algorithm to higher 
dimensions is straight forward, here we shall merely give 
an overview of the method which has been discussed in 
detail in I. 



A. QMC algorithm 

The cornerstone of the new approach is the canon- 
ical Lang-Firsov transformatioujl. which separates the 
polaron effects, due to the electron-phonon interaction, 
from the zero-point fluctuations of the harmonic oscilla- 
tors in Eq. The transformed model with one electron 
takes the form^ 

H^K + P-Ep, A- = -t^c|c^-e''^(J^'-P^) (2) 

with P as defined in Eq. , and the polaron binding en- 
ergy Ep — XW/2. The parameter 7, which corresponds 
to the displacement of the harmonic oscillator in the pres- 
ence of an electroufSi is given by 7^ = 2Ep/uj. The 
method employs a Trotter decomposition of the imagi- 
nary time axis into L intervals of size Ar = where 
(3 = (/cbT)"^ is the inverse temperature. The parti- 
tion function, obtained by integrating out the phonon 
coordinateSf2. is given by 

Zl = const. J Vp Wb Wf , (3) 

where J Dp denotes the L x dimensional integral over 
all phonon momenta p, and N is the linear size of the 
lattice in D dimensions. The bosonic weight Wb is defined 
as e"^'^'^'' with the bosonic action 

N 
i=l 
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Pi = {Pi,l^ ■ ■ ■Pi,L), and a tridiagonal L x L matrix A 
with nonzero elements 



A 



UJ 1 
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2 cjAt2 



jAt2 



(5) 



and periodic boundary conditions L+1 — 1. As discussed 
in I, the electronic weight Wi is given by 



(6) 



r=l 



Here is A' with the phonon operators pi replaced by 
the values Pi^r on the rth Trotter slice. The exponential 
of the hopping term can be written as 



DrKDl 



(7) 



{Dr)jj' ^ Sjj'e' 



where h^^ is the x A^^^ tight-binding hopping matrix 
for the lattice under consideration. In fact, this is the 
only nontrivial change compared to the one-dimensional 
case. As pointed out in I, for a single electron, the 
fermionic trace can easily be calculated from the ma- 
trix representation of fl as Wf = ^ ^ Qu . Due to the 
complex-valued hopping term, W{ is not strictly positive, 
which gives rise to the minus-sign problem discussed be- 
low. We would like to mention that, in contrast to some 
determinant QMC methods ,-34 the L-fold matrix prod- 
uct involved in the calculation of the matrix Q is well 
conditioned also for large systems at low temperatures, 
so that a time-consuming numerical stabilization is not 
necessary. 

In I, we have introduced the so-called principal com- 
ponent representation for the phonon degrees of freedom. 
In terms of the latter, the bosonic weight takes the simple 
Gaussian form 



(8) 



with the principal components = A^/'^pi. It is strictly 
positive. Finally, the efficiency of the QMC algorithm 
can be greatly improved by the use of a reweighting of 
the probability distribution. In our case, this amounts to 
transferring all the influence of the electronic degrees of 
freedom — which are treated exactly — to the observables 
calculated via 



(0) 



{Wi)h 



(9) 



with the expectation values with respect to the bosonic 
weight Wh defined as 



(0)b = 



/ Ppwb 0{p) 
J Vpwh 



(10) 



As outlined in I, the expectation values defined in 
Eq. (|10|l can be determined during the QMC simula- 
tion. The applicability of the reweighting method has 



been discussed in detail in I, and here we only mention 
that we have carefully checked that it also works reliably 
in higher dimensions. This could be expected, since the 
physics of the Holstein polaron is rather similar in all di- 
mensions (see Sec. Details about the calculation of 
observables within our formalism can be found in I. Due 
to the analytic integration over the phonon coordinates x 
used here, interesting observables such as the correlation 
functions J2i{''^i^i+s) are difficult to measure accurately. 
Other quantities such as the quasiparticle weight, and 
the closely related effective mass^^ can be determined 
from the one-electron Green function at long imaginary 
times)^ but results but would not be as accurate as ex- 
isting work on the one-electron case considered here (see, 
e.g., Refs.'28',31, andU^). Therefore, we have restricted 
ourselves to the kinetic energy of the electron, which con- 
tains a lot of information about the small-polaron cross 
over. For the more demanding many-electron case, to 
which our method can be extended (see Sec. 0), other 
methods produce far less reliable data. 

We are now in a position to discuss the actual QMC 
procedure which has been explained more thoroughly in 
I. From the above discussion and I, it is obvious that we 
only simulate the phonons, while the electronic degrees 
of freedom have been integrated out analytically. This 
is equivalent to the method proposed by Blankenbecler 
et alti^ for the many-electron case. However, here the 
Monte Carlo sampling is independent of the fermionic 
weight Wf given by Eq. lO, with the latter entering 
the calculations only via the reweighting procedure [see 
Eq. . Thereby, we avoid the problem of nonpositive 
weights during the QMC updates, since we do not use 
Wf as a probability for accepting or rejecting new phonon 
configurations. Nevertheless, the sign problem still man- 
ifests itself in terms of statistical errors, as can be seen 
from Eq. iPJ . Similar to other occurrences of the minus- 
sign problem, e.g., for the case of the Hubbard model 
away from half filling;^ simulations would become very 
difficult if {wf)h should tend to zero. 

Owing to the Gaussian form of the bosonic action, 
when written in terms of the principal components ^, the 
latter can be sampled exactly by drawing random num- 
bers from a normal distribution. In contrast to usual 
Markov chain Monte Carlo simulations, every new con- 
figuration is accepted, and no autocorrelations between 
successive values of the ^ exist. This is in strong contrast 
to conventional QMC methods for the Holstein model, 
also with many electrons (see discussion in I). Except for 
situations in which the phonons are integrated out an- 
alytically, simulations become extremely difficult at low 
temperatures and for small phonon frequencies due to 
strongly increasing autocorrelations. We regard the com- 
plete absence of such correlations as a major advantage 
of our method. After the ^ have been updated for each 
site of the L x space-time lattice, a transformation 
back to the momenta p is performed using the matrix 
Then, for each observable of interest, Owf as well 
as the fermionic weight W{ are calculated for the current 
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phonon configuration, after which the next update can 
begin. At the end of the program, results for observables 
are obtained according to Eq. © (see also I). 



B. Sign problem 

Before we come to a discussion of the sign problem 
in the approach presented here, we would like to give a 
quick review of its occurrence in other QMC methods. 
The situation is best illustrated for the case of the world- 
line algorithm (see, e.g., Ref . l37|) . The use of the latter 
to simulate systems of interacting fermions is restricted 
to one dimension by the Fermi statistics of the electrons. 
This is a consequence of the negative matrix elements w, 
which appear when two fermion world lines wind around 
each other one or more times as they traverse the space- 
time lattice. With other methods, despite the occurrence 
of the sign problem, simulations can still be carried out 
in many situations. Since the QMC sampling requires 
nonnegative weights, one uses {wl instead of w. As a 
consequence, the sign of the fermionic weight, sgn(w), is 
treated as part of the observables. The latter then have 
to be calculated via 



(C>sgn(w)}| 
(sgn(w))|^ 



(11) 



and it is obvious that simulations will become extremely 
difficult if the denominator in Eq. tends to zero. In 
fact, if a sign problem is present, the number of measure- 
ments during a QMC run has to be increased by a factor 
oc (sign)~^ to obtain results of the same accuracyi^ As 
discussed in Ref. 's^, using | w \ instead of w corresponds 
to simulating an effective model of hard-core bosons. The 
average sign can be written as 



(sign) 



-/3V(/„-/|„|) 



(12) 



where and /|^| denote the free energy per site of the 
fermionic and bosonic model, respectively, and V is the 
volume of the system. From Eq. H12|l . it is obvious that 
(sign) decreases exponentially with increasing f3 and V. 

The auxiliary-field method^ for the Hubbard model 
faces similar problems. Here the weight of a configura- 
tion is given by the product of determinants for | and 
I electrons, respectively. The product is strictly posi- 
tive only for half filling, whereas simulations for other 
particle densities become very demanding at low tem- 
peratures and/or for large systems. Since at large U 
the determinant is similar to the weight of world lines 
of Hubbard-Stratonovitch variables;^ this similarity of 
the dependence of (sign) on the parameters of the sys- 
tem is not surprising. Finally, there is no sign problem 
in determinantal grand-canonical simulations of the Hol- 
stein model at any filling, since the coupling of electrons 
and phonons is the same for both spin directions.''*' 

The sign problem in the current approach clearly has a 
fundamentally different origin, since there is only a single 



electron in the system, so that no winding of world lines 
around each other can take place. In fact, the negative 
fermionic weights are a result of the complex phase fac- 
tor in the transformed hopping term [Eq. Q]. We will 
see below that, in contrast to other methods, the sign 
problem in our approach diminishes, i.e., (sign) 1, 
with increasing system size, and a detailed investigation 
of its dependence on the parameters of the model and the 
choice of boundary conditions will be given in Sec. II VI 



C. Numerical details and performance 

As discussed in I, the error due to the Trotter decompo- 
sition is proportional to (Ar)^. We perform simulations 
at different values of At, typically At = 0.1, 0.075 and 
0.05, and exploit the linear dependence of the results on 
(Ar)^ to extrapolate to At — 0. This is a common pro- 
cedure in the context of discrete-time QMC methodsf^^ 
and allows one to remove the Trotter error if the values 
of At are sufficiently small. Moreover, a given accuracy 
can often be reached using larger values of At, compared 
to calculations which do not use the abovementioned ex- 
trapolation, thereby saving computer time. 

We conclude this section by a comparison of our ap- 
proach with other QMC methods for the Holstein po- 
laron. As mentioned in Sec. the methods of de Raedt 
and Lagendijk^^-^^-^'^ and Kornilovitc h^^i^^i^^ are based 
on an analytic integration over the phonon degrees of 
freedom. This separation of electronic and bosonic de- 
grees of freedom greatly reduces the statistical noise 
due to phonon fluctuations, which are induced by the 
electron-phonon interaction. The fluctuations increase 
noticeably with decreasing phonon frequency, decreasing 
temperature and increasing electron-phonon coupling. 
Although the Lang-Firsov transformation used here per- 
forms a very similar task, namely, to separate polaron ef- 
fects from the zero-point and thermal fluctuations of the 
free oscillators (see I), the integral over the bosonic de- 
grees of freedom is calculated with Monte Carlo, thereby 
leaving us with a residual influence of the phonons. In 
fact, the numerical effort for calculations with our ap- 
proach is proportional to the LN^^, similar to the grand- 
canonical determinant method for the Holstein model>i^ 
This could be improved to a computer time ~ N^^ by 
linearizing the exponential of the hopping term to a tridi- 
agonal matrix, which has not been done in the work pre- 
sented here. Althoug h our algorithm is not as efficient as 
the methods of Refs. 'Il'l2"l3"l4"l5"l6, in which the nu- 
merical effort is independent of N, proportional to and 
depends only linearly on D, we will see in the following 
section that it is well suited for accurate simulations on 
large clusters in one and two dimensions, and on reason- 
ably large clusters in 3D. Moreover, as we are interested 
in developing a QMC method that can also be applied 
to the many-electron system in the adiabatic regime — 
in order to study, e.g., quantum phonon effects in the 
manganites (see I) — such a decrease in performance when 
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FIG. 1: Average sign of the fermionic weight Wi as a function 
of electron-phonon couphng A in one dimension for different 
cluster sizes A'^, with parameters as indicated in the figure. 
Here and in subsequent figures, lines are guides to the eye 
only, and error bars are smaller than the symbols shown. The 
data presented in Figs. 0-0] are for At = 0.05. The inset 
shows the minimum of (sign) as a function of the system size 
N (see text). 

compared to the world-line methodaiiiiSiiii^iiSiiS, is ac- 
ceptable, since the latter cannot be applied to the many- 
electron case in more than one dimension. Addition- 
ally, the advantage of the analytical integration over the 
phonons is biggest in the one-electron case, in which they 
represent the majority of the degrees of freedom in the 
path-integral representation of the partition function^ 
This is no longer true at finite band filling, where the 
contributions of fermions and bosons are comparable. 

IV. RESULTS 

This section is divided into two parts. First, in 
Sec. IIV Al we investigate in detail the dependence of 
the aforementioned sign problem on electron-phonon cou- 
pling, dimensionality, boundary conditions, system size, 
phonon frequency, and temperature. In Sec. IIV Bl we 
present results for the electronic kinetic energy to study 
the small-polaron cross over discussed in Sec. ^ 

A. Sign problem 

Here we are interested in the average sign of the 
fermionic weight Wf, which is given by 

(sign) = J. — 1^ , (13) 

{\Wf\)h 

where the expectation value (• • •)b has been defined in 
Eq. H1U|) . In Fig. ^ we show its dependence on electron- 
phonon coupling and system size. We would like to point 
out that while the extrapolation to At = 0, discussed in 
Sec. mil has been used for the results shown in the fol- 
lowing subsection, the calculations for the average sign 
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FIG. 2: As in Fig. but for the case of two-dimensional 
clusters with (a) periodic boundary conditions and (b) open 
boundary conditions in real space. The insets show the min- 
imum of (sign) as a function of the linear system size N. 



have been performed for a single, fixed value At = 0.05, 
for which the Trotter error is smaller than statistical er- 
rors. Some other results for (sign) in one dimension have 
already been reported in I. 

From the general discussion in Sec. IIIII it is clear that 
the sign problem encountered in the present approach is 
of a different nature than in, for example, QMC simula- 
tions of the Hubbard model. As reported in I, for the Hol- 
stcin polaron problem under consideration, it is most pro- 
nounced for small systems, low temperatures and small 
phonon frequencies u) <C 1 . Therefore, the bulk of results 
presented below will be for such a set of "worst case" 
parameters, including N = A, f3t = 10, and ui — 0.1. Fig- 
ure^ shows that, in one dimension, the average sign of 
Wf in the critical region of intermediate electron-phonon 
coupling increases quickly as the system size increases 
from = 4 to = 16, which is in strong contrast to 
Eq. H12|l. This increase of the minimum as a function of 
A^ is also shown in the inset of Fig. As we have only 
calculated (sign) for a finite number of A-values, an ap- 
proximation for the minimum has been determined using 
a spline interpolation. 

The minimum of (sign) occurs near A = 1, where the 
transition from a large to a small polaron takes place (see 
Sec.nil. The increase of statistical errors in this regime 
as a consequence of the sign problem is similar to the situ- 
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FIG. 3: Dependence of the average sign of W{ on (a) phonon 
frequency ui and (b) inverse temperature /3 on a 4 x 4 cluster. 
The insets show the minimum of (sign) as a function of (a) 
phonon frequency Co, and (b) inverse temperature /3t. 




FIG. 4: Dependence of the average sign of Wf on the hnear 
system size A*' in three dimensions. The inset shows the min- 
imum of (sign) as a function of the linear system size A'' for 
two different temperatures. 



ation encountered with simulations of the untransformed 
model (see I). However, the use of the transformed model 
still gives significantly more accurate results for the same 
number of measurements, in particular for low tempera- 
tures and small phonon frequencies. Finally, one may be 
tempted to explain the unusual system-size dependence 
of the sign problem by ascribing its origin to the periodic 



boundary conditions in real space. If the latter were in- 
deed the source of the sign problem, the boundary effects 
would decrease with increasing system size, in accordance 
with the results of Fig. This possibility has been in- 
vestigated, and we shall see below that the sign problem 
persists also for the case of open boundary conditions. 

We now come to the Holstein model in two dimensions. 
In Fig.|2Ia), we show results for (sign) for different lattice 
sizes, again starting with a very small linear dimension 
N. All other parameters are the same as before, in par- 
ticular pt — 10 and u! = 0.1. Obviously, for the smallest 
cluster size shown, the minimum of the average sign has 
diminished to a value of approximately 0.1, so that large 
numbers of measurements are necessary. However, simi- 
lar to one dimension, (sign) increases with increasing sys- 
tem size, and for the largest system size shown (N — 12), 
we find a rather uncritical minimum value of about 0.5. 

The results for open boundary conditions, shown in 
Fig-EIb), reveal that for small clusters, the average sign 
increases compared to the case of periodic boundary con- 
ditions. However, with increasing system size, (sign) 
quickly converges to the same values, independent of the 
boundary conditions. This is just what one would ex- 
pect, since with increasing N, the effect of the choice of 
boundary conditions on the properties of the system di- 
minishes. Moreover, we can conclude that the negative 
weights do not simply result from hopping processes of 
the electron across the periodic boundaries, since in that 
case we would expect (sign) = 1 for open boundary con- 
ditions, in contrast to Fig. |2Ib). Similar behavior has 
been found in one and three dimensions. 

The influence of the phonon frequency to on the average 
sign is shown in Fig. I^a) . Clearly, the sign problem is 
most noticeable for small values of Q, while it diminishes 
quickly as we increase the adiabatic ratio (see inset). This 
is very similar to the small-polaron cross over. As dis- 
cussed in Sec. IHl the latter sharpens significantly with 
decreasing to, while an abrupt transition is completely 
absent in the nonadiabatic regime lu > 1. Moreover, the 
coupling Ac at which the minimum of (sign) occurs in- 
creases with Q, in agreement with the aforementioned 
small-polaron condition AD/iD > 1. 

In Fig. I^Ib), we report the average sign as a function 
of A, and for different inverse temperatures p. Again we 
have taken iV = 4 and cZ; = 0.1, the parameters for which 
the sign problem is most noticeable. While for (3t = 10, 
the minimum of (sign) lies below 0.1, the situation is 
much better already for (3t — 5, sls shown in the inset. 
At even higher temperature /3t = 1, the fermionic weight 
is always positive so that we have (sign) = 1 for all A. 
The dependence of the sign problem on temperature is 
therefore similar to other QMC methods [see Eq. (|12|) ]. 
although we do not find a simple exponential relation. 

Finally, we also present results for (sign) in three di- 
mensions, for lattices of different linear size TV, and again 
for the parameters Cu = 0.1 and f3t = 10 as a function of 
A. Figure 01 reveals that the minimum of (sign) in 3D 
has an even more pronounced form than in two dimen- 
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sions. The sign problem diminishes shghtly as we increase 
the system size from = 4 to = 6. However, accu- 
rate simulations in this regime are still quite demanding. 
As we will see in Sec. IIVBI a temperature oi pt — 5 
is sufficient to study the low-temperature properties of 
the Holstein model. For this higher temperature, and all 
other parameters unchanged, (sign) is close to 1 even for 
TV = 4 (see inset in Fig. , so that accurate results can 
be obtained even for small phonon frequencies in three 
dimensions (see Sec. 11 V B)! . 

In summary, the investigation of the sign problem has 
shown that our method works well for a large range of val- 
ues of phonon frequency, electron-phonon coupling and 
temperature, as long as the system size is large enough. 
This contrasts, for example, with the world-line algo- 
rithm of Kornilovitchji^ii& which is restricted to interme- 
diate and strong coupling, as well as low temperatures 
and rather large phonon frequencies (see also Sec. ITHl . 
also by a minus-sign problem. Although we have inves- 
tigated the influence of all important parameters on the 
sign problem, a physical interpretation of its origin has 
not emerged. Nevertheless, it is clear that the negative 
fermionic weights are a result of the phase factors in the 
Lang-Firsov transformed hopping term in Eq. 10) . This 
is similar to the sign problem which occurs, e.g., in sim- 
ulations of electron-phonon models in an external mag- 
netic field, but here the phonon fields pij vary with 
time (l) and position (i), and are coupled in imaginary 
time by the bosonic action [see Eqs. Q and Q]. More- 
over, the dependence of the sign problem on uj, A, /3 and 
N bears striking resemblance to the influence of these 
parameters on the properties of the Holstein polaron. 
In particular, its reduction with increasing system size 
may be a consequence of the dilution of the system (the 
particle density n ^ as A'^ — > cxd, in the one-electron 
case) . Finally, it is interesting to note that the large sta- 
tistical fluctuations — resulting from the sign problem in 
the case of the transformed model — occur at exactly the 
same points in parameter space as in the untransformed 
model. This suggests a strong correlation between the 
minimum in (sign) and the transition to a small polaron, 
which both occur at or near Ac , similar to simulations of 
other models, in which the sign problem occurs exactly 
where the most interesting physics is going on, i.e., in the 
vicinity of phase transitionsi^ 
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FIG. 5: Normalized kinetic energy E-^ [see Eq. 1141 1 as a 
function of electron-phonon coupling A for different values of 
the phonon frequency on a 12 x 12 lattice. The results 
shown in Figs. IKl- llll have been obtained by extrapolating the 
QMC data to At = (see text). 

with i^k = 1 for T = and A = 0. Due to the 
large amount of work that has been devoted to the Hol- 
stein polaron in the past (see Sec. ^ lot is known 
about the transition from a large to a small polaron 
with increasing electron-phonon coupling. In this work, 
we therefore concentrate on those aspects which have 
not been studied in a systematic way so far. To this 
end, we exploit the advantages of our QMC method 
which allows us to investigate, in particular, finite-size 
and finite-temperature effects. The latter have only 
been touched upon briefly by previous authors employing 
Q]y[C,iiJ^ii3Jiii5J4 who focused On very largeii^'^^Oi 
or infinite systems in the ground statei^^ or at 
two different temperatures^i*i^ However, since a large 
amount of work has been done using ED or other meth- 
ods based on diagonalization of small clusters, it is es- 
sential to study the convergence of the results with in- 
creasing system size. Additionally, we shall also present a 
comparison of the kinetic energy in one, two and three di- 
mensions. As pointed out before, the results shown here 
have been obtained by extrapolating to At = 0, thereby 
eliminating the error due to the Trotter discretization. 

1. Two dimensions 



B. Holstein polaron 

The most interesting observable, which is easily acces- 
sible with our method (see discussion in Sec. IHf and 
allows us to investigate the small-polaron cross over, is 
the one-electron kinetic energy i?k ~ (K), given by the 
expectation value of the first term in Hamiltonian 
In order to compare results for different dimensions, we 
define the normalized quantity 

= ^^k/(-2iD) (14) 



To study the small-polaron cross over, previous authors 
focused on the effective mass of the electroni^iSSiSSi^ and 
the quasiparticle weight i^^iSSi^ Unfortunately, as pointed 
out in Sec. 1111 Al these observables cannot be calculated 
directly with our QMC algorithm. Nevertheless, we can 
examine many interesting aspects of the cross over by 
calculating the kinetic energy of the electron, and com- 
parison will be made to existing workiiitiSiiiiiS^ 

We begin with the dependence of the cross over on 
the phonon frequency. To this end, we present in Fig. [S] 
results for the kinetic energy calculated for N = 12, 
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FIG. 6; Normalized kinetic energy Ey^ as a function of 
electron-phonon coupling A for different linear dimensions A'' 
of the system. 



f3t — 10 and different values of lu. The large range of 
the adiabatic ratio in Fig. [S] 0.1 < w < 8.0, shows the 
ability of our method to give accurate results for almost 
arbitrary values of the phonon frequency, especially in 
the adiabatic regime (D 1. This is in contrast to the 
method of Kornilovitch;i^ii& which is restricted to (D > 1 
and A > 1 by a severe minus-sign problem. In our case, 
the only limitations regarding the accessible values of uj 
are the moderate sign problem, discussed in Sec. IIV Al 
which for small systems and low temperatures gives rise 
to a noticeable increase of statistical errors as u) — > 0, and 
the increasing Trotter error as (D — > oo, which requires the 
use of more and more time slices (see discussion in 1). 

Figure m also shows the well known fact that the tran- 
sition from large to small polaron near A = 1 sharpens 
considerably with decreasing phonon frequency. While 
there is an abrupt decrease in E'k for a) = 0.1, the cross 
over is very smooth for u) > 1. Additionally, we see from 
Fig. El that the cross over position Ac increases with tu in 
the nonadiabatic regime. The physics of the transition 
to a small polaron has been discussed, e.g., by Capone et 
al^ In the adiabatic regime (w <C 1), the cross over is en- 
tirely determined by the balance of kinetic and electron- 
phonon coupling- or displacement energy, where the lat- 
ter is given by the polaron binding energy Ep. As soon 



FIG. 7: Normalized kinetic energy as a function of the 
inverse of the linear size A'^ of the system, and for different 
values of the electron-phonon coupling A. As indicated in 
the legend, some curves have been shifted, in order to allow 
for a better representation. All curves are monotonic within 
statistical errors. 



as the gain in displacement energy outweighs the loss in 
kinetic energy, the electron localizes in a potential well 
and forms a polaron. The parameter A is defined as the 
ratio of these two contributions and may be written as 
A = Ep/{W/2) {—W/2 is the kinetic energy of a free 
electron at T = 0). Therefore, in the adiabatic regime, 
the cross over occurs at Ac = 1. With increasing w, the 
lattice energy becomes more and more important, since 
more energy is required to excite phonons. As a conse- 
quence, the distortions of the lattice around the position 
of the electron — giving rise to the large effective mass 
and low mobility in the small polaron regime — are much 
smaller, and the local oscillators will predominantly be 
in their ground state. Thus, even for A > 1, where a 
trapped state is energetically favored in the adiabatic 
regime, the electron remains mobile. The decrease of i?k 
with increasing A is a result of the exponentially decreas- 
ing overlap of a displaced and an undisplaced harmonic 
oscillator in its ground state, which reduces the hop- 
ping matrix element between neighboring lattice sites. 
In the nonadiabatic regimeja small polaron is formed if 
7^ — Ep/uj > 1 (see Ref. l33|) . where 7 is the induced 
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FIG. 8: Normalized kinetic energy E]^ as a function of 
electron-phonon coupling A for different inverse temperatures 
/3. 



lattice distortion at the site of the electron which follows 
from the Lang-Firsov transformation.^ This is identical 
to the condition AD/iD > 1 given above. The larger lat- 
tice energy also gives rise to the more gradual decrease of 
E]f^ for intermediate and large values of Co. In particular, 
the kinetic energy is much larger for A > 1 and lu > 1 
than for a; < 1. 

Finally, like in one dimension,^ Ey^ remains small but 
finite even at very strong coupling, which is a conse- 
quence of the fast and undirected motion of the electron 
inside the polaron. As pointed out by Kornilovitch^ 
the kinetic energy in the small polaron regime is there- 
fore not related to the effective mass of the electron, the 
latter exhibitin g an exponential decrease as a function of 
A (Refs.lHllili and ,31). 

To address the issue of finite-size effects, we have calcu- 
lated i^k for pt = 10, ct) = 0.1 and different linear lattice 
sizes N = 4-12 [see Fig. Ela)]. The choice of (I> = 0.1 
is reasonable since the large polaron, which exists for 
A < Ac, is most extended for small phonon frequencies, 
as discussed in Sec. so that finite-size effects can be 
expected to be largest. To illustrate this point, we also 
present results for a larger phonon frequency Cu — 2.0 
[Fig-EIb)]. In the latter case, the local oscillators can re- 
spond very quickly to the motion of the electron, and the 
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FIG. 9: Normalized kinetic energy E]^ as a function of 
temperature and for different values of the electron- 

phonon coupling A. As indicated in the legend, some curves 
have been shifted, in order to allow for a better representa- 
tion. 



extension of the phonon cloud or lattice distortion sur- 
rounding the electron is much smaller. Since for lu > 1 
the transition to a small polaron happens at larger A, the 
polaron will be larger for intermediate values 1 ^ A < 2 
compared to the adiabatic regime, in which Ac = 1. In 
total, we therefore expect to have smaller finite-size ef- 
fects for uj — 2.0 than for a) = 0.1 as long as we are in the 
weak-coupling regime A < 1 , while the opposite should be 
true for 1 < A < 2. For A > 2, a highly immobile polaron 
state exists in both cases and results should therefore be 
virtually independent of N. All this is well confirmed by 
the results shown in Fig. |^ 

Figure|n|also reveals that the results begin to saturate 
for > 8, as pointed out previously by KornilovitchiM 
However, in contrast to the one-dimensional casCf^ where 
we performed simulations for A^ as large as 32, we find 
a nonnegligible dependence on A'^ up to the largest sys- 
tem size (A^ — 12). This is better illustrated in Fig. [7| 
in which we show Ey^ as a function of 1/N, and for sev- 
eral values of A. To allow for a better representation, 
some curves have been shifted, as indicated in the leg- 
end. From Fig. [7| we see that Ey changes very little for 
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N > A for the case of strong coupling A = 2 and, in fact, 
remains constant within the error bars, similar to the ID 
results in I. For smaller values of A, no such saturation 
is found on the scale of Fig. [7| and the behavior of the 
kinetic energy for large N is almost linear when plotted 
as a function of 1/A^. We used a linear fit of the data 
for = 8, 10, and 12 to obtain an approximation to the 
thermodynamic limit. The results are shown in Fig. 
Obviously, Ey^ for iV = oo has decreased noticeably for 
small values of A (including A = 0), while it remains 
almost unchanged in the small polaron regime. The de- 
crease of i5k for A < Ac can easily be understood if we 
consider the fact that our method works at a finite tem- 
perature 1/(3. As a consequence, for very small N, the 
energy gap between the ground state with fc = and the 
first excited state with fc 7^ is larger than the thermal 
energy (/3t)~^. With increasing system size, thermal pop- 
ulation of excited states becomes possible. For TV large 
enough {N sa 20 for /?t = 10 and A = 0), results converge 
to those for N — (x and, in fact, the extrapolated data 
for A = 0, shown in Fig. agree well with the results 
for a free electron on an infinite lattice. We ascribe the 
smallness of this finite-temperature effect in the strong- 
coupling regime above Ac to the very small width of the 
polaron band (see, e.g., Ref. I2T). Consequently, the low- 
energy coherent states with different k have very similar 
energies. 

Figure Eta) also shows that for small phonon frequen- 
cies, the influence of the lattice size is largest near Ac. 
The cross over is sharper for small systems than for the 
case of large N. A similar behavior has been found by 
Marsiglio for the one-dimensional model. '^^ 

Finally, as pointed out by Kornilovitch,!^ despite the 
good convergence of the results for quantities such as the 
kinetic energy even for small ~ 4, other observables, 
such as the effective mass, rely on the knowledge of the 
energy of states with infinitely small momenta. There- 
fore, they will show a more pronounced dependence on 
N when calculated on small clusters accessible with, e.g., 
ED. 

We also investigated the effect of temperature on our 
results, again for a) = 0.1 and ui = 2.0. The results, 
shown in Fig. |H| indicate that E]^ is more affected by 
the finite temperature of the simulation in the adiabatic 
case oj = 0.1 [Fig.|H|[a)]. This is a consequence of the fact 
that calculations at finite temperatures only give ground- 
state-like results when (3oj 3> 1. Clearly, this condition 
is much more difficult to meet for uj = 0.1, and requires 
larger values (3t > 10. 

The changes of -Ek with temperature result from an 
interplay of several effects. For A = 0, the kinetic en- 
ergy approaches its full noninteracting value of —2tD 
(i.e., E'k = 1) as T ^ 0. At finite T, however, states 
with nonzero total quasimomentum k will contribute and 
thereby lead to a decrease of E]^. As discussed above, this 
effect of temperature on Ey^ is expected to decrease with 
increasing A, and to be extremely small in the strong- 
coupling regime. In the adiabatic case [Fig. |HIa)], the 



transition at Ac = 1 is smeared out at high tempera- 
tures. For both uj = 0.1 and oj = 2.0, a qualitative 
change in behavior occurs near Ac- Ey decreases with 
increasing temperature for A < Ac, whereas the oppo- 
site is true for A > Ac. The behavior above Ac can be 
understood by considering the electronic hopping ampli- 
tude, given by the overlap of the wave functions of a 
displaced and an undisplaced oscillator at neighboring 
sites. While the latter is exponentially reduced with in- 
creasing electron-phonon coupling at T = (see above), 
it increases with temperature since the oscillators can 
occupy excited states, corresponding to wave functions 
which are more spread out than the ground state. More- 
over, the thermal energy allows the electron to overcome 
the potential barrier more easily. This thermally acti- 
vated hopping — which shows up more clearly for small 
phonon frequencies since, in this case, phonon excitations 
require less energy — has been studied already a long time 
ago by Holstein?^ 

It is interesting to notice that the dependence of Ey 
on temperature, as shown in Fig. O is almost linear at 
low temperatures and for uj = 2.0 [see Fig.lSJ^b)], while in 
the adiabatic case [Fig. Eta)] this is only true for strong 
electron-phonon coupling. However, as pointed out be- 
fore, for w = 0.1, an inverse temperature f3t — 10 is not 
sufficient to obtain well-converged results. Therefore, a 
linear dependence, as for uj — 2.0, may still be found at 
lower temperatures. Such calculations, with the accu- 
racy of the results presented here, would be more time- 
consuming. The data for N = 12, P = 10, At = 0.05, 
and A = 1.0, i.e., in the cross over regime where statistical 
errors are largest (about 0.5%), required about 28 days of 
CPU time on an Intel Xeon 2600 MHz computer. How- 
ever, as has been mentioned previously, the numerical ef- 
fort may be reduced by a factor of order by linearizing 
the exponential of the hopping matrix, denoted as k (see 
Sec. nil CI ) Similar to the finite-size scaling performed 
above, we also linearly extrapolated the data for pt = 10 
and 15 to the zero-temperature limit pt = 00. The re- 
sults for oj = 2.0, for which such a linear scaling is rea- 
sonable, are shown in Fig.Ob). While the general trend 
agrees well with our expectations based on the finite- 
temperature data shown, the scaling procedure clearly 
overestimates the temperature effects, thereby leading to 
spurious values -Ek > 1 at A = 0. However, this can easily 
be understood keeping in mind that results saturate at 
low-enough values of pt, so that the linear extrapolation 
used here becomes insufficient. 

Finally, we would like to mention that the results of de 
Raedt and Lagendijl4ii*iSii^ were also given for pt = 5 
(and UJ — 1), but the small number of Trotter slices {L — 
32 in their case) gives rise to relatively large systematic 
errors,^^ This is not the case for the results of Ref. O in 
which the same extrapolation to At = was employed 
as here. Finally, the method of KornilovitcbMii^ is free 
of such errors, but only permits one to calculate ground- 
state properties for a restricted range of Co and A. 
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FIG. 10: Normalized kinetic energy i5k as a function of 
electron-phonon coupling A for different linear dimensions A*' 
of the lattice. 
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2. Three dimensions 

In contrast to the two-dimensional case dis- 
cussed above, less work has been done in three 
dimensions liiiiSiiiiSiiSiSSi^SiSi In fact, we are only aware 
of one calculation of the kinetic energy, which is by de 
Racdt and Lagendijk (DRL)J^ To compare with their 
work, we chose the same values for the phonon frequency 
w = 1.0 and temperature pt = 5. As pointed out in 
Sec. mil the numerical effort for calculations with our 
method, which is proportional to N^^ for the algorithm 
in the form used here but could be reduced to iV^^, 
restricts us to smaller systems than those considered by 
DRLJiiSii^ For simplicity, we have therefore limited 
ourselves to a maximum of = 6, for which results can 
easily be obtained within a reasonable amount of com- 
puter time, while the data presented in Refs. Illll2ll3l 
is for N = 32. To be more specific, our calculations for 
one value of A, for iV = 6 and Ar = 0.05, took about 
10 h on an Intel Xeon 2600 MHz computer. Due to the 



relatively small system size in our work, it is important 
to study to what extent the results are converged with 
respect to N. To this end, in Fig. ^1 we present i?k as 
a function of A for iV = 4, 5, and 6. Surprisingly, the 
results are already satisfactorily converged, keeping in 
mind the rather small values of N. There is a maximal 
change of less than 20% in the transition region at A = 1, 
while i?k remains almost constant for small and large 
A, as the linear size increases from = 4 to = 6. 
Thus, increasing N further will not change the results 
qualitatively, although the finite temperature of our 
simulations will manifest itself in a way similar to the 
two-dimensional case. Our findings agree well with the 
results of DRL.^^ The main difference is that for weak 
coupling, our results are closer to the zero-temperature 
values (e.g., i^k = 1 at A = 0). The reason for this 
discrepancy — despite the fact that we have used the 
same temperature — is the smaller lattice size in our 
calculations. In contrast to the two-dimensional case 
considered above, we have not performed a scaling to 
A^ = oo in 3D, since the clusters under consideration are 
too small to reveal a systematic power law dependence 
as a function of 1/A^. 

Finally, wc wish to investigate the effect of dimension- 
ality on the small-polaron cross over. Therefore, we com- 
pare i?k in one, two and three dimensions using f3t — 5, 
A^ = 6, and uj = 1.0. The dependence on D, which is 
shown in Fig. 1111 is in perfect agreement with previous 
work. The transition from a mobile large polaron to a 
small polaron, moving in a very narrow but still coherent 
band, sharpens considerably with increasing dimension 
of the system, and while Ey^ only displays a gradual de- 
crease in ID — without any signs of an abrupt change at 
A — 1 — we find a sharp and well-defined transition in 
three dimensions. 

We conclude this section by comparing the accuracy of 
our results with the QMC methods of DRLfiiii2iiii4 and 
Kornilovitchii^ii& As discussed in Sec. lIIII their main ad- 
vantage over our method is the fact that they allow one to 
obtain data which are essentially free of finite-size effects, 
in any dimension D = 1-3, and with modest computa- 
tional effort. However, we have seen above that even in 
three dimensions, where the limitation of our algorithm is 
most noticeable, results are reasonably converged. While 
the approach of Refs. and is limited to T = 0, 
DRL's method as well as the current approach can, in 
principle, be used to study any temperature. Apart from 
the sign problem discussed in Sec. 1111 lj| the only limi- 
tation which occurs is the fact that one has to increase 
the number of Trotter slices as /3 ^ oo, so as to keep 
the Trotter error smaller than the statistical errors. This 
situation can be greatly improved by extrapolating re- 
sults to Ar = 0, as has been done in this work, while 
the continuous-time algorithm of KornilovitchiSii& is free 
of any such discretization errors. The accuracy of the 
results presented here depends on Co, (3t, N and A. Away 
from A ~ 1, error bars are use usually smaller than the 
line width, corresponding to relative errors of less than 
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0.5%. This is comparable to the accuracy of the resuhs 
given by KornilovitcbtS and significantly more accurate 
than the original results by DRLiii*i^ 

V. CONCLUSIONS 

Extending the work of Ref. to the case of the two- 
and three dimensional Holstein model with one electron, 
we have shown that our QMC approach allows accurate 
calculations with modest computational effort for a large 
range of parameters. In particular, the minus-sign prob- 
lem, due to the Lang-Firsov transformation, has been 
found to diminish quickly with increasing system size, 
and not to have a significant effect on simulations. We 
have presented a detailed study of the small-polaron cross 
over and its dependence on the parameters of the sys- 
tem. In particular, we have focused on finite-size and 
finite-temperature effects, which have not been investi- 
gated systematically before. 

As discussed above, our approach is not as fast as other 
QMC methods for the Holstein polaron,iUaiiiiliSii& the 
main limitation being the restriction to smaller but still 
reasonably large lattices. This difference in performance 
is acceptable keeping in mind that it can be extended 
to the many-electron case (see below), in contrast to 
the world-line methods of Refs. illiil^l3il4il5iil6» . which 
work best for the case of a single electron or two elec- 
trons of opposite spin, and face a severe sign problem 
in more than one dimension, similar to other world-line 
methods^ Finally, despite the sign problem, we can per- 



form accurate simulations for virtually all physically in- 
teresting values of ui and A, in contrast to the method of 
Refs. Eland EE 

Motivated by the promising results of this work and 
Ref. 0, our next objective will be a generalization to the 
many-electron case. To this end, it is important to no- 
tice the striking similarity of the QMC method presented 
here with the determinant method first introduced by 
Blankenbecler et al.^ which has been successfully used 
to study superconductivity and charge-density-wave for- 
mation in the Holstein model liSiiSiSSiSiiS^ In fact, the 
main difference between our one-electron algorithm and 
the corresponding grand-canonical method is the calcu- 
lation of the fermionic weight, which in case of the latter 
is given by the determinant of l + fl. Here 1 denotes the 
unity operator, and is the corresponding generalization 
of Eq. © to many electrons. Consequently, the numeri- 
cal effort is expected to be similar to the simulations with 
one electron. However, an important open question con- 
cerns the dependence of the sign problem on the number 
of electrons. This work is currently in progress. 
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